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ABSTRACT 

We describe implementation and tests of sink particle algorithms in the Eu- 
lerian grid-based code Athena. Introduction of sink particles enables long-term 
evolution of systems in which localized collapse occurs, and it is impractical (or 
unnecessary) to resolve the accretion shocks at the centers of collapsing regions. 
We discuss similarities and differences of our methods compared to other imple- 
mentations of sink particles. Our criteria for sink creation are motivated by the 
properties of the Larson-Penston collapse solution. We use standard particle- 
mesh methods to compute particle and gas gravity together. Accretion of mass 
and momenta onto sinks is computed using fluxes returned by the Riemann solver. 
A series of tests based on previous analytic and numerical collapse solutions is 
used to validate our method and implementation. We demonstrate use of our 
code for applications with a simulation of planar converging supersonic turbu- 
lent flow, in which multiple cores form and collapse to create sinks; these sinks 
continue to interact and accrete from their surroundings over several Myr. 

Subject headings: Hydrodynamics — (magnetohydrodynamics:) MHD — meth- 
ods: numerical — ISM: clouds — stars: formation 



1. Introduction 

Gravitational collapse is a common feature of many gaseous astrophysical systems, and 
specialized methods are required in order to follow collapse in time-dependent hydrodynamic 
simulations. These numerical issues are particularly important in studies of star formation. 
As gravitational collapse develops, material converges to a central point from all directions 
to create a large density peak. For simulations that follow the formation of multiple self- 
gravitating prestellar cores within large-scale clouds, the true structures that form as a 
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consequence of core collapse are generally so small that the profile surrounding each collapse 
center becomes un-resolvable for grid-based codes, even if adaptive mesh refinement (AMR) 
is adopted. For example, stellar radii are ~ 10 11 cm, whereas that of a giant molecular 
cloud (GMC) is ~ 10 20 cm. With a dynamic range > 10 9 , it is not possible to spatially 
resolve a central post-shock protostar at the same time as capturing the large scale flows 
that lead to its formation, when multiple collapse centers are simultaneously present. When 
the central density in a collapsing region becomes too large, gradients in gas pressure and 
gravity from the central cell to the neighboring cells cannot be resolved, such that correct 
mass and momentum fluxes cannot be computed by the numerical solvers. A simulation 
cannot continue under these conditions. 

The huge dynamic range involved in gravitational co llapse can also lead to diffi culties 
due to the time step restriction by the Courant condition (IRichtmyer &: Mortonlll994l ). The 
time-scale set by self-gravity varies oc p _1//2 , so that an increase by a factor of > 10 6 in 
density relative to ambient conditions (as would apply within the centers of collapsed cores) 
implies a decrease in the time step by a factor ~ 10 3 . 

A practical way to deal with the above difficulties, for numerical models that are fo- 
cused on large-scale dynamics, is to establish a minimum spatial resolution and introduce 
sink particles. When gravitational collapse occurs, the unresolved high density peaks are 
eliminated from the grid and replaced with sink particles. After a sink is created, the mate- 
rial can flow smoothly toward the center of collapse, with the profile remaining well resolved 
near the sink. Subsequent to the creation of sinks, gas and sink particles are integrated 
simultaneously, including mutual gravitational forces. Provided that the flow onto sinks is 
supersonic, introducing them will not affect the dynamics of the upstream flow. Accretion 
onto sinks should be implemented in a way that conserves total mass and momentum of the 
system. 



Bate et al.l (119951 ) firs t introduced sink parti cle t echniques in a smoothe d particle hy- 



drodynamics (SPH) code. iKrumholz et al.l ( 120041 ) and iFederrath et al.l ( 120101 ) implemented 
sink particles in the grid-based codes ORION and FLASH, respectively, and extensively 
discuss tests of their methods. In the past several years, similar im plementations hay e 
been made for a number of other well-established c odes, such as ENZ O rtWang et al.lboiol) , 



RAMSES (Dubois eta 



Krumholz et al.l ( 120041 1. 



2010; iTevssier et al.l 120111). and GAD GET (Ijappsen et al. 



Wang et al.l ( 12010luDubois et al.l ( l2010luTevssier et al.l(l2011f l adopt 



2005 



the same methodology, in cluding the criteria for cr e ation of sink particles and the accretion 

1 ) use d an early version of the im- 
[ 2010 ) . In t h ese implementations. 



rate onto sink particles. IVazquez-Semadeni et al.l (120 



plementation in FLAS H described by 



Federrath et al 



Krumholz et~all ( l2004h . IFederrath et al.l ( hoioh . IWang et all feoioh . IPadoan fc Nordlund 



-3 - 



( 2011 ) 



Vazguez-Semadeni et al.l ( 120111 ) mainly focus on star formation simulations; iDubois et al. 



Teyssier et al.l (j201l[ ) create sink particles to replace super-massive black holes in cos- 



mological simulations. 

In this paper, we pr esent details of our implementation and tests of s i nk particles in the 
grid- based code Athena (IGardiner k Stonell2005l . 120081 ; IStone et al.ll2008l ; IStone k Gardiner 
20091 ). In Section 2, we begin by introducing the Eulerian code and outlining our methods for 
implementing sink particles. We physically motivate criteria we adopt for creating sinks, and 
describe our methods for treating gravity, accretion, merging, and orbit integration of sink 
particles. In Section 3, we present a series of tests of our methods. These include orbits of 
two particles, collapse of self-gravitating spheres with a range of initial conditions (including 
self-similar solutions, and converging supersonic flows). We demonstrate Galilean invariance 
of our methods. To illustrate the capabilities of our methods for typical applications, in 
Section 4 we consider evolution of a turbulent medium with a large-scale supersonic flow 
converging to create a dense slab. We follow the fragmentation of the slab into multiple 
cores, and the subsequent evolution of the system as sink particles are created and grow. 
Finally, Section 5 summarizes our presentation. 



2. Numerical Methodology 

2.1. Athena Code 

For the simulations presented in this paper we use the three-dimensional ( 3D) code 
Athena JGardiner k Stone! I2OO5I . I2OO8I : Istone et alibood : Istone k Gardiner! |2009[ ). Athena 
is a grid-based code that uses higher order Godunov methods to evolve the time-dependent 
equations of compressible hydrodynamics and magnetodynamics (MHD), allowing for self- 
gravity, ra diative heating and cooli ng, and other microphysics, on either a Cartesian or 
cylindrical ( jSkinner k Ostrikerll2010l ) grid. In this paper, we only refer to the hydrodynamics 
and self-gravity capabilities of the code, with Cartesian coordinates. The hydrodynamics 
equations solved are the mass, momentum, and energy equations: 



dp ^ 
— + V 

dt 



d(pv 

dE 
~dt 



and the Poisson equation, 



(pv) = 0, 

+ V-{ P VV + P*) = -V(<K + $ ex t), 

+ V- [(E + P)v) =pv V($ + $ cxt ), 
V 2 $ = AnGp, 



(1) 

(2) 
(3) 

(4) 
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where P* is a diagonal tensor with P* = P I , and P is the gas pressure, P is the total energy 
density 

E = - At + V, (5) 
7 — 1 2 

$ is the gravitational potential of the gas, and <3> ext is an external gravitational potential. In 
this paper, we shall consider isothermal flows, in which the energy equation ([3]) is replaced 
by the relation P = (? s p. Here cj. = kT/fi is the square of the isothermal sound speed, for T 
the temperature and \i the mean mass per-particle. 

To solve the Poisson equation under varying boundary conditions, we adopt two different 
gravity solvers. For periodic boundary conditions in all directions, we use the fast Fourier 
transformation (FFT) method. For boundary conditions that are perio dic in-plane (x — y) 



and o pen in the ^-direction, we use the FFT method developed by iKoyama &: Ostriker 



(120091 ). For open boundaries in al l direc tions, we use a new solver based on the method 



described in iHockney &: Eastwood! (Il98ll ). which employs FFTs in a zero-padded domain 



eight times as large as the computational box (see Appendix). 

Athena also includes static mesh refinement (SMR), but we do not refer to these capa- 
bilities in the current work. Implementation of our sink particle algorithm with SMR will 
be discussed in a future publication. 



2.2. Creation of Sink Particles 



Different criteria for the creation of sink particles have been discussed in the past decade . 
A density threshold (IJappsen et al.l 120051 ; iPadoan fc Nordlundl l201ll ; iKrumholz et al.l 120041 ) 
is the simplest criterion for creation of a sink particle. High density regions may form in as- 
tronomical systems as a consequence of different physical processes, most importantly strong 
supersonic shocks and gravitational collapse. If a high density region is not gravitationally 
bound, it might subsequently be destroyed by large scale motions that induce rarefactions. 
Thus, care must be taken to select an appropriate density thre shold, and to include a dditional 
criteria that must be satisfied before creating sink particles ( IFederrath et al.ll2010l ). 



Using AMR simulations, iTruelove et al.l (119971 ) showed that grid-scale numerical noise 
could grow to cause artificial self-gravitating fragmentation if the local Jeans scale (Lj = 
c s {tx /Gp) 1 / 2 ) is not resolved by at least four cells. This criterion gives a density threshold 



7T 



PTr 



16 GAx 2 



(6) 



for Ax the simulation cell size. 
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The implementation oflKrumholz et al.l (120041) se ts the threshold density fo r sink p article 
creation in a cell to p^. iBanerjee et al.l (120091 ) and I Vazquez- Semadeni et al.l (120 111 ) adopt 
criteria that first checks whether density exceeds a threshold den sity and th e n che cks if 
the zone in que s tion is a local gravitational potential minimum. iBate et al.l (119951 ) and 



Federrath et al.l (120101 ) perform a series of checks including a density threshold check, a 
converging flow check, a local gravitational potential minimum check, and additional checks 
that evaluate whether a region is strongly self-gravitating. The converging-flow check of 
Federrath et al.l (120101 ) requires that the flow surrounding a candidate zone is converging 
along all directions, stricter than the condition V • v < 0. 

For simulations involving gravitational collapse, a sink particle should only be cre- 
ated at the center of a region that is collapsing. As the sink particle creation criteria of 
Krumholz et al.l (120041 ) do not specifically limit particle creation to a single collapse cen- 
ter, they adopt the approach of merging the spurious sink particles within a given collaps- 
ing region so that one final sink particle is created inside each potential well. The extra 
cri teria — checks f o r a lo c al potential minimum o r a g ravitationally bound state - — ado pted 
by IBanerjee et al.l (120091 ) , Federrath et al.l ( 120101 ) and IVazquez-Semadeni et al.l (120111 ) fur- 
ther limit the initial creation of sink particles, so that only one sink particle is c r eated for 
each local potential minimum. We shall adopt similar criteria to Federrath et al.l ( 120101 ) to 
ensure that single sink particles are created inside regions that are collapsing. 



2.2.1. Density threshold 



Our choice of density threshold, p t h r , for sink particle creat ion is motivate d by the well- 



known solution for self-gravitating collapse first obtained by iLarsonl (119691 ) and iPenston 
(119691 ) (hereafter LP). For collapse of an initially-static, gravitationally-unstable isothermal 
sphere, LP showed that a singular density profile 



P LP ( r ) 



8.86c;: 
4ttGV 2 



(7) 



is reached. Numerical studies with a range of initial conditions have shown that the LP 



asymptotic solution is in fact an "attractor" for isotherma 



core collapse, no matter how the 



collapse is initiated 


Bodenheimer & Sweieart 


1968: Larson 


1969 


; Penston 


1969; 


Hunter 


1977; 


Foster & Chevalier 


199: 


5: 


Oeino et al. 


199C 


); 


Hennebelle et al. 


2003 


; iMotovama & Yoshida 


2003; 


Vorobvov & Basu 


2005; 


Gomez et al. 


2007; 


Burkert & Alves 


2009 


) . Some of the above 



models start from static unstable configurations, and others from static, stable configurations 
that are subjected to an imposed compression, either from enhanced external pressure or a 
converging velocity field, or a core-core collision. These models all show that collapse starts 
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from outside and propagates in, and that the central velocity is comparable to the value 
— 3.28c s at the time of singularity for mation, when the density prof ile approaches the inverse- 
square LP asymptotic solution p L p. Gong fc Ostriker (j2009 . 2011 ) showed that the collapse 
of cores forming inside converging supersonic flows also approaches the LP solution, whether 
the flow is sp herically symmetric or is turbulent, with no special symmetry. iGong fc Ostriker 
( 120091 . 120111 ) also showed that the duration of collapse (starting when the core is ~ 10 times 
the ambient density) is typically a free-fall time at the mean core density. 

When collapse occurs, the LP profile would result in a density 



p LP (0.5Ax) 



8.86 



7T 



GAx 2 



at a distance Ax/2 from the center. We take this as the density threshold (pthr) for sink 
particle creation. We note that this value is 14.4 times the value pir of Equation ([6]). A 
potential concern is that this might lead to artificial fragmentation. However, we find (see 
section 4.3) that using the threshold from Equation ([8]), no additional particles are created 
compared to cases in which we instead adopt Equation <Q (pn-) for the sink particle density 
threshold. After a sink particle is created, all cells centered at r < 2Ax become part of a 
"sink region" (see below), such that the density of cells exterior to the sink region satisfy 
the Truelove criterion, with p L p(2Ax) = 8.86c^/(167rGAx 2 ) < px r - 

Our standard choice of density threshold is given by Equation ([8]). However, we have 
also tested other choices, as discussed in Section 4.3, and found similar results. 



2.2.2. Control volume 

Surrounding each sink particle is a sink control volume where the gas flow cannot be 
resolved. As a sink particle moves from one cell to another, the control volume moves with 
the sink particle, such that the sink particle is always located within the central zone of 
the cubic control volume. The sink control volume is generally set to (3As) 3 although we 
find similar results for value (5Ax) 3 and (7Ax) 3 . The effective radius of the control volume 
is ~ r ctrh with r ctr i = 1.5 Ax for a control volume (3Ax) 3 . As gravity is unresolved at 
the same scale, the same control volume is adopted for the particle-mesh algorithm we use 
to compute the sink particle's gravity (see section 2.4). Once a control volume has been 
created, it acts similar to ghost zones bounding the simulation domain. At every time step, 
the density, momentum, and energy of the cells inside each control volume are reset using 
outflow boundary conditions from the active grid (i.e. via extrapolation from surrounding 
non-sink zones). 
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2.2.3. Gravitational potential minimum check 



If the density of a cell with integer indexes (i,j,k) exceeds p tnr , and its distance to 
any existing sink particle is larger than 2r ctr i, a "temporary" control volume is created 
surrounding it. This temporary control volume has same size as the real control volume for 
a sink p article. We then check if the central cell is a potential minimum inside this control 
volume (IFederrath et al.ll2010l ). If the potential minimum test is satisfied, we apply further 
tests. 



2.2.4- Converging flow check 



As in IFederrath et al.l ( 120101 ) but less restrictive, we test whether the candidate sink 
cell at the center of the control volume is surrounded by a converging flow: V ■ v < 0. 
Under most circumstances, the tests for high density, a gravitational potential minimum, 
and a converging flow would guarantee the region surrounding this cell is under gravitational 
collapse. However, under the special circumstance of a strong shock produced by a converging 
flow, all of the above criteria might be met, but a region would still disperse i f the converging 
flow were not sustained until the region becomes gravitationally bound ( IFederrath et al. 
20101 ); this would occur when the volume of the high-density post-shock region becomes 
sufficiently large. 



2.2.5. Gravitationally bound state check 

As a last check before sink particl e creation, we test whe ther the total energy inside the 
temporary control volume is negative (IFederrath et al.ll2010f ): 

-Egrav + E th + E kin < 0, (9) 

where E grav is the gravitational potential energy, E t h is the thermal energy, and -Ekin is the 
kinetic energy. 

The control-volume gravitational energy is calculated as 

£grav = ^2p(i,j,k)A$(i,j,k), (10) 

ijk 

where A$(i,j,k) = $>(i,j,k) — $0 is the potential difference between the cell with integer 
index (i,j,k) and the potential at the "edge" of the local potential well, $o- For $ , we 
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compute the average value of the potential in all zones immediately outside the temporary 
control volume. 

The thermal energy and kinetic energy are calculated as follows: 

E t h = -^2p(i,j,k)c 2 s (i,j,k), (11) 

ijk 

E km = ~^2p(i,j,k)\v(i,j,k) - v cia \ 2 ; (12) 

ijk 

here v cm is the velocity of the center of mass of the control volume. 

If a cell passes all the checks above, a sink particle is created and a permanent control 
volume is tagged around it. The initial mass and momentum of the sink particle are set by 
the sums over all zones within the control volume. Note that the cells within the sink control 
volume are not modified either at the moment of sink creation or by subsequent accretion 
because they are effectively ghost zones; their values are updated at the same time as other 
boundary conditions (see Section 2.2.2). Also, as noted above, the sink control volume is 
redefined whenever the sink particle moves to a new zone within the computational grid. 
However, the particle can move within a given cell without redefining the sink control volume. 



2.3. Gas Accretion Onto Sinks 



A key aspect of any si nk particle implementa tion is to ensure that the gas accretion 
rate onto sinks is accurate. iKrumholz et al.l (120041 ) use the Bondi-Hoyle accretion formula 
to approximate the accretion rate, with the sound speed and flow velocity set by host cell 
values. The density in the Bondi-Hoyl e formula is set based on the mean density in a 
local accretion zone (see their section 2.4). iFederrath et al.l (120101 ) take a simpler approach, 
removing mass from any cell within the control volume to the sink particles if the density in 
that cell exceeds pxr- 

Motivated by the concept of the sink region as "internal" ghost zones, in our algorithm 
the accretion rates of mass and momentum to each sink particle are calculated based on 
the fluxes returned by the Riemann solver at the interfaces between the control volume and 
the surrounding computational grid. Tests to confirm that the accretion rate agrees with 
analytic solutions are discussed in Section 3. As sink particles cross the border of one cell 
to enter the next cell, the mass and momentum differences between the new and old ghost 
zones are combined with the fluxes returned by the Riemann solver to conserve the mass 
and momentum of the whole simulation domain, including both gas and particles. 
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2.4. Integration and Merging of Sink Particles 



Once they have been created, sink particle posi tion and veloci ties must be integrated in 
time. To do this, we use the leapfrog method (e.g. 
At given by: 

At\ 



Springell 120051 ). with updates over time 



or 



U(At) = D 



U(At) = K 



2 j K(At)D 



At 



A A 

"2"; 



(13) 



(14) 



Here D(At) and K(At) are the "drift" and "kick" operators respectively, and D(At) updates 
a particle's position without changing its momentum, while K(At) does the opposite; U(At) 
is the time evolution operator for an interval At. Both drift-kick-drift (DKD, Equation f[T3T) ) 
and kick-drift-kick (KDK, Equation f[T4|) ) schemes are implemented. For all the simulations 
presented in this work , the KDK sch eme is adopted because it is superior to the DKD scheme 
for variable time step (jSpringelll2005l ). The time step of the whole simulation is also restricted 
by the sink particle velocities: sink particles cannot travel further than one grid zone in one 
time step. 

Velocity updates of the sink particles are set based on the gravitational field at the 
particle's (smoothed) location. This gravitational field must include a contribution from 
all the o ther particles, as well as fro m the gas. We use the triangular-shaped-cloud (TSC) 
scheme (IHockney &: Eastwood! Il98ll ) to calculate the gravity produced by each particle as 
well as the force each particle feels. In this method the mass of each sink particle is smoothed 
to the n^ gc cells surrounding it, where n TSC is the number of cells used in one direction. The 
weights along the x direction are expressed as follows: 



W?{Ah) 



2(n„ 



2(" a 



-I) 2 



-2#) 



( n TSC 


-1)2 


2(n Tgc — | 


o Ah \ 
Z Ax> 


(n Tgc - 


1)2 , 







l TSC 1 I 1 
2 ' 



•,-1 



: 



(15) 



/ — 1 n TSC 1 _ 1 

) 1 ■*•)■••> 2 



n, 



l TSC ^ 



( n TSC 

(1+2H) 2 /_ 
2(n TSC -l)2- 

Here Ah is the x-distance between the sink particle and the center of the cell where it resides 
along the x dimension. The index I ranges from — (n TSC — l)/2 to {n TSC — l)/2 and the cell 
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with I = is where the s ink particle resides. Taking n Tgc = 3, Equation f|T5|) reduces the 



TSC weights described in iHockney &: Eastwood! (Il98ll ): 



2Ah\ 2 3 {Ah\ 2 „„ 1 / 2A/A 2 



= o i - -it— ; w<f = - - hr- ; w+i = - l + — — . (16) 

Analogous weights are defined along the y and z- dimensions, giving and W*. 

The weight of the mass in any of the n^ gc zones of the smoothing volume surrounding 
the particle is a product of the weights from the x, y, and z dimensions. After the weights 
of a sink particle are computed, the particle's mass m p is applied to the grid. Within each 
particle's smoothing volume, the particle's effective density at zone (l,m,n) is given by 
W \ W xAyAT V ■ Outside of the sink's control volume, the density is equal to that of the gas. 
After the combined gas+particle density is defined, the solution $ of the Poisson equation 
is obtained via FFTs (using the FFTW package, see www.fftw.org). For each zone within 
the sink control volume, the gravitational field / lmn is computed taking differences of the 
potential. The gravitational force on each sink particle is computed using the same weights 
as in Equation (TT5]) : 

f= E E E W l XW Xfl™, (17) 
lmn 

where l,m and n run from — (n TSC — l)/2 to (n TSC — l)/2. 

With a TSC smoothing volume of n^, gc zones, the effective radius for smoothing by the 
TSC algorithm is r TSC = n c Ax/2. For nxsc = 3,5 and 7, r TSC is 1.5Ax, 2.5Ax and 3.5Ax 
respectively. With larger r TSC , sink particles' masses are smoothed to a larger volume, and 
the gravity is softened more. A smaller value of r TSC gives more accurate gravity near a 
particle. However, the hydrodynamic fluxes are better resolved if the radius r ctr i of the sink 
control volume (see section 2.2.2) is larger. We have tested different combinations of r ctr i 
and r TSC with r ctr i > r TSC for the test problems described in Section 3.1 and 3.2. We find 
that different r ctr i give nearly the same accretion rate. To maximize resolution, we therefore 
adopt n TSC = 3 and r ctr i = r TSC = 1.5Ax as our standard choice. 

Sink particles in our algorithm represent unresolved star/disk/envelope structures. Be- 
cause the gas structure and gravity is not resolved within the control volume of each sink 
particle, the detailed interaction of star/disk/envelope systems that collide or pass near 
each other cannot directly be followed. We therefore take a conservative approach and 
merge two sink particles whenever their control volumes overlap. A new sink particle is 
created at the center of mass of the two previous sink particles, with total mass and momen- 
tum conserved. We note that provided the effective mass distribution of each sink remains 
spherically-symmetric and concentrated at a very small scale, it would be safe to adopt 
less-strict merger criteria in methods that compute interparticle forces directly rather than 
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using particle-me sh gravity; in this case close particle-particle interactions require sub-time- 
stepping (see e.g. iKrumholz et al.l 12004 ; iFederrath et al.ll2010f ). Because torques on the gas 
are not well resolved when the sink region has only 3 3 zones, we do not presently track the 
spin angular momentum of the sink particles. 



3. Tests of the Method 

3.1. Particle Orbits 

To test the TSC algorithm and the leapfrog integrator, we consider a problem in which 
two particles with equal mass orbit their common center in circular orbits. The initial circular 
speed of the two particles is v = ^^2Gm/ d, where G is the gravitational constant, m is 
the mass of each sink particle, and d is the distance between two sink particles. The gravity 
between sink particles and gas is disabled, so that sink particles only feel gravity from each 
other. Figure [1] shows the trajectory of one particle for 10 orbits for d/L = 0.2 (top panels) 
and d/L = 0.3 (bottom panels) with resolutions 32 3 , 64 3 and 128 3 from left to right. Here, L 
is the simulation box size, corresponding to either 5 or 3.3 times the interparticle separation 
d (which is the semimajor axis for the reduced-particle Kepler orbit). Vacuum boundary 
conditions are used for the gravity solver. 

As long as the gravity is well resolved, the combination of the TSC scheme and leapfrog 
integration gives accurate orbits. For d/L = 0.2, the radius of the orbit is 3.2Ax, QAAx and 
12.8Aa: for resolution 32 3 , 64 3 and 128 3 respectively. For 32 3 resolution with r ctr i = 1.5Ax, 
the cubic smoothing volumes for each particle are separated by only 2 — 3 zones, such that 
the gravitational potential is not well approximated by a point mass. For resolution 64 3 and 
128 3 , this distance ranges (8 — 10)Ax and (21 — 22)Ax respectively. We find that for 64 3 , 
the orbits are a bit coarse, whereas for 128 3 , the orbits are nearly perfect. Correspondingly, 
for the d/L = 0.3 cases, the orbit radius is 4.8Ax,9.6Ax and 19.2Ax at resolution 32 3 ,64 3 
and 128 3 , with circular orbits well represented at the higher resolutions. We conclude that 
provided the orbit radius is resolved by > 10 zones, quite accurate orbits are obtained, with 
somewhat lower quality orbits at smaller separations. However, even for orbit radii as small 
as 3 zones (or ~ 0.03pc for typical simulations of core/star for mation in molecular clo uds), 



orbits are approximately circular. By comparison, we note that (IFederrath et al.ll201ll ) have 
found that resolution of a vortex in their hydrodynamic simulations on a Cartesian grid 
requires a radius of 15 grid cells. 
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3.2. Self-similar Collapse of Isothermal Spheres 



To confirm that our algorithms yield the correct accretion rate, we compare our numer- 
ical results with analytic solutio ns. In part icular, we compare to the family of self-similar 
accretions solutions analyzed by IShul ( 119771 ). The density profile for a generalized singular 
isothermal sphere is 

Ar 2 

*> = (18 > 

For A > 2, equilibrium is not possible because gravity exceeds the gas pressure gradient 
everywhere; for these initial conditions, a sphere would globally collapse everywhere. The 
case A = 2 with zero velocity everywhere corresponds to an unstable hydrostatic equilibrium. 



Shul (119771 ) analyzes a family of self-similar solutions in which the density profile in the outer 
parts approaches Equation (TTSl . while the inner part approaches a free-fall profile with 



p oc r 



" 3 / 2 . The outer part has v oc — r _1 , while in the inner part the velocity approaches 
free-fall v oc — r _1//2 . For any value of A, the central accretion rate is constant, such that the 
central mass M oc t, and we can define a dimensionless accretion rate M/(c 3 /G). 

We h ave tested a series of values of A ranging from 2.0004 to 4, the same as the values in 
Table 1 of Hfi93k0 The initial density and velocity profiles are the solution of Equations 
(11)-(12) in Shu ( 19771 ). obtained using a four-step Runge-Kutta integration scheme. For 
code units, we adopt an arbitrary density po, together with length scale, Lj = c s (tt / 'Gpo) 1 ^ 2 
and time scale tj = Lj/c s . For all of the tests in this section, the simulation domain size is 
(4Lj) 3 , and the resolution is 129 3 . Vacuum boundaries are used for the gravity solver. 

The initial radial density and velocity profiles are truncated at r max = l.BLj. Outside 
this radius, the initial density is set to p(r = 1.5Lj), and the initial velocity is set to 
zero. To convert the self-similar solutions to initial density and velocity profiles input to the 
simulation, we choose an initial time t/tj = 0.43 such that in the case A = 2.0004, the initial 
radius of the expansion wave is 11% of the box size and 29% of the initial radius of the sphere. 
In the case A = 2.0004, the total mass of the sphere within r = 1.5Lj is 18.9c 3 (47rG 3 po) -1 ^ 2 - 
A sink particle is introduced at the center of the sphere at the beginning of the simulation, 
with it s ini t ial m ass set to the mass of the initial profile within r = r ctr i (Equations (8), 
(10) in IShul (119771 )). Note that the uniform density outside the sphere leads to "outside- in" 
collapse because of the imbalance of gravity and gas pressure at the truncation radius. This 
process will not affect the initial "inside-out" accretion to the central sink particle, however. 



Krumholz et al.1 (|2004l ) show (their Fig 1) a compariso n of the numerically-computed den sity and velocity 



profile s to the semi-analytic collapse solution of lShul (|19770 for the singular isothermal sphere: iFederrath et al 



(|2010h show (their Fig. 10) the density and velocity profiles of their numerical solution at several times for 
the case A=29.3. 
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Figure [2] shows the accretion history of the central sink particle for the cases A = 2.0004 
(crosses) and 4 (open circles). The solid lines are the analytic solutions from Shu (1977). 
Before effects from the outer boundary conditions begin to affect the accreting region, the 
simulation evidently reproduces the analytical solutions extremely well. For the case A = 
2.0004, the particle mass at t — 6(47rGpo) -1//2 > near the end of the linear stage, has reached 
46% the sphere's initial mass, increasing by a factor three from its initial value. For each 
value of A, we compute the accretion rate in the simulation using a fit to the sink particle's 
mass versus time during the stage when this evolution is linear. 

Figure [3] shows the accretion rates for models based on Shu' s solutions with different 



values of A. The line shows the values from the analytic solution of IShul (119771 ). and the solid 
black dots are based on measurements from our numerical simulations. The 3D simulations 
reproduce nearly exactly the predicted accretion rate from the analytic solutions. 

Figure H] shows a sample cross-section of the density and velocity field near the center 
of the collapsing singular sphere for A = 2.0004. Since A is close to 2, the simulated sphere 
recovers Shu's famous "expansion-wave" solution: the outer part retains a static singular 
isothermal equilibrium profile before the expansion wave arrives, and the inner part is free- 
falling towards the center. Figure [5] shows the radial density profiles and velocity profiles 
from the simulation at different instants during the collapse. The filled circles are the profiles 
from the numerical simulation and the solid lines are the corresponding analytic solutions. 
Our simulation reproduces the analytic solutions. The "expansion wave" is clearly seen from 
both the density and velocity profiles. In both density and velocity profiles, it is seen that 
gas very near the boundaries collapses, and the outer density is slightly altered from p oc r~ 2 . 
This behavior is due to conditions near the boundaries, but this process does not affect the 
collapse in the interior. 



3.3. Galilean Invariance of Accretion 

To confirm that the accretion is properly computed for moving particles, we consider 
tests in which both the particle and surrounding gas are initialized with a bulk flow across 
the grid. The initial conditions of these models are exactly same as the A = 2.0004 case in 
Section 3.2, but with an additional uniform bulk flow of speed v^d = 0, 0.5c s , 1.5c s and 2.5c s . 
If the update of the mass and momentum of the sink particle is correct, it will continue to 
move with the surrounding gas sphere at the same velocity. The mass accretion rate onto 
the sink particle should also be the same for different bulk flow speeds. For this problem, 
periodic boundaries are adopted for both gas and the gravity solver. For all of the tests in 
this section, the simulation size is (4Lj) 3 , and the resolution is 129 3 . 
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We note that the boundary conditions in Section 3.2 are open, different from what we 
adopt in this section. Here, we adopt periodic boundary conditions so that the sphere may 
pass through the boundary of the simulation domain (for large Mach number cases). Because 
of the difference in boundary conditions, the accretion rate is different from that in Section 
3.2; however we shall show that the accretion rates are all consistent with each other for 
different Mach numbers. 

Figure E] shows the speed and mass of the sink particle versus time, for all tests. For the 
case of zero bulk advection, the speed of the particle remains zero all the time. For u a( j ^ 0, 
the mean speed of the sink particle agrees with the bulk speed, with small oscillations. 
During an interval Ax/vad, the time for the sink particle to cross a zone, the control volume 
does not move. However, as discussed in Section 2.2.2, the sink control volume must be 
reset when the particle crosses a zone boundary. During the time interval Ax/v^, the 
mass and momentum fluxes into the control volume are not exactly symmetric, because 
the gravitational potential is not symmetric if the particle is offset from its zone center. 
The downwind material enters the control volume with relatively larger momentum than 
the upwind material, and addition of net positive momentum accelerates the sink particle. 
The speed of the sink particle therefore temporarily increases. As the particle crosses the 
border of one cell to enter the next cell, the sink control volume is shifted, and mass and 
momentum differences from cells entering and leaving the control volume are applied to the 
particle to conserve mass and momentum. As a consequence, the particle's speed is reduced 
immediately after crossing a cell edge. At late times, the fractional change in the momentum 
becomes small, so the oscillation amplitude drops. 

Figure [6] shows that the mass of the sink particle increases at the same rate for all the 
tests, at varying t> ad . At early stages, before material originating near the sphere's boundaries 
accretes to the sink particle, differences in the accretion rate are completely negligible, and 
and differences remain very small even after material originating near the boundaries reaches 
the center. 

We conclude that our sink particle algorithm satisfies Galilean invariance, with the 
accretion rate the same whether or not the flow and the grid are in relative motion. 



3.4. Collapse of a Bonnor-Ebert Sphere 



To test our algorithm for the the creation of a sink particle, together with the a ccretion 
after format i on, w e run 3D simulations of the evolution of a Bonnor-Ebert sphere (IBonnor 
19561 ; lEbertl Il955l . hereafter BE) and compare the results to a ID spherically-symmetric 
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simulation conducted with another code. We consider a BE sphere with a radius slightly 
larger than the critical radius, r = l.lr B E,crit for r B E,crit = 0.274c s /(7r/Gp e ) 1 / 2 , where p e is the 
edge density. The density inside the sphere is everywhere twice as large as the equilibrium 
value. Exterior the sphere, the density is set to 0.001p e . The initial velocity is everywhere set 
to zero for the initial conditions. The gas boundary conditions are outflow. The resolution 
is 129 3 , and the box size is [0.301c s / (n/Gpe) 1 ^} 3 . 

Figure [7] shows the cross-section of the density and velocity field across the center of the 
BE sphere immediately before the central sink particle is created. Evidently, the inner part 
collapses inwards, following the "outside-in" collapse described by |L arson! ( 119691 ). Because 
the outer boundary of the sphere is not in equilibrium (for this isothermal simulation, we 
do not have a hot confining medium) , the outer part of the sphere expands outwards at the 
same time. Only the very late accretion to the sink particle is affected by the early expansion 
of the outer parts of the sphere. 

Figure |8] shows the accretion rate and the mass of the sink p a rticle from our 3D simu- 
lation, in comparison with the ID simulation of iGong fc Ostrikerl ( 120091 ) obtained with the 
ZEUS code. The solid lines are from the 3D simulation and the dashed lines are from the 
ID simulation. The accretion rates immediately after the creation of the particle differ, but 
the accretion rates at later stages are almost exactly the same. The peak accretion rate is 
higher for ID than 3D because it is measured closer to the sink particle (0.5 zones vs. 1.5 
zones). In addition, the point of singularity formation cannot be resolved as well on a 3D 
Cartesian grid as ID spherically-symmetric grid. However, it is evident that the 3D code 
captures the overall accretion history very well. 

To further check the results of the model evolution, Figure [9] compares the density and 
velocity profiles for 3D and ID simulations at the instant of the sink particle creation. The 
density profiles both approach the Larson- Penston singular solution (Equation ([7])), while 
the velocity profiles approach the limit v = — 3.4c s , in the inner part. 

Figure [10] shows the radial density and velocity evolution for the 3D simulation for this 
test. The solid lines show the profiles before the sink particle is created, the dashed lines 
show the profiles at the instant of the sink particle creation, and th e dotted lines show the 
profiles during the envelope infall stage (see IGong fc Ostrikerl 120091 ) . Profiles are separated 
by time At = 0.043(7r/Gp e ) 1 / 2 . The collapse of the inner part (and the expansion of the 
outer part) is clearly seen. The density approaches the singular LP profile p oc r~ 2 (Equation 
((7])) at the instant of collapse, and the corresponding velocity profile is flat in the inner part, 
approaching — 3.4c s . After the creation of the sink particle, the inner density and velocity 
profiles approach free-fall. 
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3.5. Converging Supersonic Flows 



Gong fc Ostrikerl (120091 ) presented a unified model for dense molecular core formation 
and evolution, based on spherically-symmetric simulations. In that work, dense cores are 
not present as either stable or unstable density concentration in the initial conditions, but 
are built by the convergence of supersonic flows. Post-shock compressed gas accumulates 
over time in stagnant, shock-bounded regions. When the core accumulates enough mass, 
it becomes gravitationally supercritical and collapses, leading to formation of a protostar 
at the center. This is followed by a stage of infall of the envelope onto the protostar, and 
subsequent accretion of ambient material. 

To ensure the accretion rate of sink particles is accurate for cores formed by supersonic 
flows, we run a simulation of a 3D spherical convergin g supersonic flow, and compare to 



the accretion rate obtained by iGong fc Ostrikerl (120091 ) using a ID spherically-symmetric 



simulation. The initial density is uniform everywhere, with a value po (this represents a 
typical density within GMCs). The flow in the initial conditions converges to the center 
everywhere at Mach 2. The size of the simulation box is 1.6c s / (AttG po) 1 ^ 2 . This size is 
chosen so that it is large enough for the post-shock compressed region to grow until it 
collapses. The simulation is run with 129 3 cells. For the ID spherical symmetric simulation, 
the initial density and velocity profiles are the same as 3D but the model is run with 64 
zones. For the ID simulation, a sink cell is introduced at the center after collapse occurs, 



with an outflow boundary condition next to the sink cell (IGong fc Ostrikerl 120091 ). 



Figure [TT] shows the accretion rate and the mass of the protostar versus time during 
the envelope infall stage, comparing the 3D and ID simulation results. The accretion rate is 
nearly same for the two cases, and the mass history of the "protostar" is comparable. For the 
ID simulation, the accretion rate will decrease to exactly the inflow value 87rGpo r box at late 
time. For 3D simulation, however, we cannot create an ideal spherical converging flow given 
our cubic domain, such that the accretion rate decreases after t ~ 0.2('47rG ? / 3n)~ 1 ^ 2 - Note that 



the pe ak of the accretion rate for the ID model is smaller than the value in lGong fc Ostriker 



(120091 ) for Mach 2, due to lower resolution here. 



4. Planar Converging Supersonic Flow with Sink Particles 



The most generic configuration for converging supersonic flow is planar, and lGong fc Ostriker 



(120 111 ) presented a set of 3D numerical simulations with this geometry - com bined with turbu- 
l ent p erturbations - to study the core building and collapse. The models of IGong fc Ostriker 
(120 111 ) represent a localized region within a giant molecular cloud, in which there is an over- 
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all convergence in the velocity field (produced by the largest-scale motions in the cloud), 
combined with a turbulent power spectrum i n which linewidth incre ases with size. Since 
sink particle techniques were not employed in iGong fc Ostrikerl (1201 if ) , each model simula- 
tion was stopped when the most evolved core collapsed. To demonstrate the capabilities for 
following multiple core collapse an d evolution when s ink p articles are introduced, here we 
rerun a sample simulation from the IGong fc Ostrikerl (120111 ) suite. 



As in lGong fc Ostrikerl ( 120111 ). we adopt the isothermal approximation. The isothermal 
sound speed at a temperature T is 



0.20 km s 



-i 



T 



10 K 



1/2 



(19) 



If the density within clouds were uniform, the spatial scale relevant for gravitational insta- 
bility would be the Jeans length 



7T 



Gp 



1/2 



2.76 pc 



10 2 cm 



-1/2 



T 



10 K 



1/2 



(20) 



evaluated at the mean density po- The corresponding Jeans mass is 



Mj = p L 3 j = c 



71 ' 



1/2 



72 M m 



™h,o 



10 2 cm 



-1/2 



T 



10 K 



3/2 



(21) 



The Jeans time at the mean cloud density is 



tr= — 



7[ 



1/2 



3.27t ff (p ) = 1.4 x 10 7 yr 



™H,0 



10 2 cm 



-1/2 



(22) 



Since we consider a planar converging flow, another relevant quantity is the total surface 
density integrated along the z direction, 



p(x, y, z)dz 



for S ee p Lj = 9.49 M ^{T /\mfl 2 (n Hfi /W cui 



p dz 
Po Lj 



(23) 



-3U/2 



For this test, the supersonic flow converges to the central plane from +z and —z di- 
rections at Mach number Ai = 5, for total relative Mach number 10. For both the whole 
domain initially and the inflowing gas subsequently, we apply perturbations following a 
Gaussian random distribution, with a Fourier power spectrum of the form 



(\6v k \ 2 ) oc k~ 



(24) 
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for \kL/2ir\ < N/2, where N is the resolution and L is the size of the simulation box 
in x an d y. The power spectru m is appropriate for supersonic turbulence as observed in 
GMCs (IMcKee &: Ostrikerl 120071 ) . The perturbation velocity fields are pre-generated with 
resolution 256 3 in a box of size Lj. The perturbation fields are advected inward from the 
^-boundaries at inflow speed M.c s : at time intervals At = Az/(A4c s ), slices of the pre- 
generated perturbation fields for v x , v y and v z are read in to update values in the ghost zones 
at the z-boundaries. 



We set the amplitude of the turbulent power spec trum at the large s cale t o Svm(L j) = 
1.3c s , which corresponds to the low amplitude case in iGong &: Ostrikerl ( 1201 ll ). The resolu- 



tion is N x x N y x N z — 256 x 256 x 96, with domain size L x x L y x L z /L 3 j = 1 x 1 x 0.375. 



4.1. Structure Evolution 

Figure [T2] shows evolution of the surface density (Equation ( 12"B"|) ) projected in the z- 
direction after the most evolved core collapses and creates a sink particle (marked as "1" in 
the images). The instants for the four images from top left to the bottom rights are: 0.301tj, 
0.349tj, 0.398tj and 0.446t j. The time interval between these images is 



At = 0.048tj = 6.72 x 10 5 yr ( ™ H '° - ) ' . (25) 

\10^ cm^/ 

The black dots and numbers mark the sink particles formed prior to the time of each snapshot 
for the first three plots. In the last frame, the magenta curves show the trajectories of sink 
particles and the black triangles show where these sinks were created. Over time, some of the 
sinks merge, and the large white solid dots show the final set of post-merger sink locations 
at t — 0. 446to- The numbers marking sinks in the images indicates the sequence of their 
creation. In the top right figure, the sink number 7 is missing since it has merged with the 
sink 1. Similarly, in the lower left, several of the sinks have already undergone mergers, and 
their corresponding numbers are not shown (5, 6, 7, 10, 11). 

Filamentary features dominate the moderate-density structure in all images. These 
structures grow from the initial turbu l ent pe rturbations, which are then amplified by self- 



-1/2 



gravity (see Fig.l in IGong &: Ostrikerl (120111 )). The localized collapse of these filamentary 



structures leads to the formation of protostars. Filaments also become more stratified over 
time as they acquire more material and contract perpendicularly under their self-gravity. 

The bottom right panel in Figure [T2] shows the trajectories of all sinks, as well as their 
merging history. There were 12 sinks created during this simulation, and four survive at 
final time. All of the other 8 have merged with other nearby sinks due to close encounters. 
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4.2. Sink Particle Mass Evolution 

Figure [13] shows the mass of sinks versus time. The solid lines show the four sinks 
that survive up until late stages. The dotted lines are sinks that undergo mergers with 
other more massive sinks. Each sink forms at the center of a dense core, and the solid 
dots show the gravitationally bou nd core masses calculated using the GRID-core finding 



algorithm ( iGong fc Ostrikerl 120111 ) immediately before the formation of the corresponding 
sink particle. The mass of every sink particle grows smoothly to reach and exceed the bound 
core mass, and keeps increasing until it merges with other sinks. As indicated in the figure, 
the final masses of sink particles are all much larger than the bound core masses at the initial 
instant of sink formation, which are in the range m core = 0.02 — O.lMj. 

The formation of sinks is divided into three groups. Protostars 1, 2, 3, 4 and 5 form at 
the earliest time; sinks 6, 7, 8, 9 and 10 form a bit later; and sinks 11 and 12 form during 
the late accretion stage. Sinks forming at earlier stages are more likely to survive to the 
end, and their masses grow significantly via mergers and late accretion. For example, the 
merging history for protostar lis: 7 — >• 1 ■<— 4 <r- 11, the merging history for protostar 2 is: 
10 — 2 ^ — 12, and the merging history for protostar 3 is: 6 — > 5 — > 3 <r- 8. 

At late stages, the masses of these surviving sink particles are very high: M\ = l.llMj, 
M 2 = 0.79Mj, M 3 = 1.13Mj and M 9 = 0.38Mj. These correspond to tens of M if 
n n,o — 10 2 — 10 3 cm -3 . In part, these sinks may end up with very high mass because 
the turbulence has low amplitude and is purely decaying. As consequence, matter is not 
prevented from accreting into the potential wells that develop. In reality, outflows, radiative 
feedback, and other energy injection would limit gas accretion onto protostars. In addition, 
some stars could be ejected due to close gravitational encounters, before acquiring a high 
mass. These physical issues will be addressed in a future publication; it is straightforward 
to implement localized feedback with rates set by the mass, age, accretion rate, etc, of each 
sink. Here, our goal is simply to test the proper implementation of sink particle algorithms 
and to demonstrate that these enable robust long-term evolution. 

We note that for the present algorithm (and similarly for other sink particle implemen- 
tations), the minimum particle mass depends on the grid resolution, on the minimum density 
threshold for sink creation, and on the size of the sink particle control volume. The sink 
particle density threshold we adopt here is p t h = P LP (Q-5Ax) = 8.86c 2 / (nG(Ax) 2 ). With 
resolution Ax = Lj/N, p t h = 8-86po{N/7r) 2 , and the mass within the central cell of the sink 
region is 0.89Mj/iV = iV _1 65 M (n H ,o/lO 2 )" 1/2 (T/lOK) 3 / 2 . The total initial sink particle 
mass are larger, due to non-negligible density in the surrounding sink cells. We find (see Fig- 
ure [12]) that the initial sink masses in this test simulation are in the range 0.02M/ — O.lMj. 
With N = 256, this ranges from ~ 6 — 30 times the minimum mass. 
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4.3. Criteria for Sink Particle Creation 

As discussed in Section 2.2.1, different density thresholds have been adopted for sink 
particle creation by different groups, motivated by different physical and computational 
considerations. Figure [TH shows a comparison of sink particle mass versus time for different 
sink particle creation criteria. The solid lines are based on criteria of exceeding the LP density 
threshold (Equation (151)) and satisfying the local potential minimum check. The tracks 
marked by pluses adopt additional criteria: the converging flow check and gravitationally 
bound state check. The dotted lines adopt the Truelove density instead of the (higher) LP 
density as a threshold, and apply a check for a local gravitational potential minimum. We 
find no differences between the sink masses for the first two sets of criteria. For the third 
set of criteria, the sink masses at birth are much lower than in the first case, because of the 
lower density threshold. Compared to initial sink masses of ~ 0.02M/ using the LP density 
criterion, initial sink masses are ~ 0.005Mj using the Truelove density criterion. However, 
these sinks evolve to follow tracks identical to those found with the other two sink criteria 
choices. Analogous comparisons at different Mach numbers show that no additional particles 
are created by artificial fragmentation when using our standard criteria, even though the LP 
density threshold is higher than the Truelove density threshold. 

We conclude that a proper density threshold and local potential minimum check are 
probably sufficient criteria for sink particle creation. By choosing a density threshold given 
by the LP solution, and requiring that the zone is at a local potential minimum, we can 
ensure that the region surrounding a cell is under gravitational collapse. 



We note that in the study of [Truelove et al.l (119971 ). there was no implementation of 



sink particles. There, the authors showed that artificial fragmentation may occur in AMR 
simulations if the ratio of the Jeans length to the grid scale is too small. However, in simu- 
lations where sink particles are introduced, the formal requirements on res olution need not 



be identical to the case whe n there are no sink particles. Previous authors (IKrumholz et al. 



2004 ; iFederrath et al.ll2010l ) showed that they were able to avoid artificial fragmentation by 
capping the density at the Truelove value (see Equation E]) • However, we find that when a 
slightly higher density criterion (motivated by LP collapse) is adopted for sink particle cre- 
ation, there is still no artificial fragmentation, and in fact the mass evolution follows similar 
tracks as it does when the Truelove density threshold is used (see Fig. IT4"|) . The similarity 
among the mass evolution tracks for different sink particle creation criteria suggests that the 
class of sink particle methods currently in use provides reliable results for core collapse and 
accretion. 



5. Summary 



We have presented an implementation of sink particles to the grid-based Eulerian hy- 
drodynamics code Athena. A standard particle-mesh method is adopted to calculate gravity 
forces by and on the sink particles, with the Poisson equation solved via FFT methods. 
We use the mass and momentum fluxes from the Riemann solver to update the mass and 
momentum of sink particles. Criteria for sink particle creation are similar to those used 
by other authors, although we suggest that a higher density threshold is motivated by the 
Larson-Penston profile that is known to develop as a generic stage of self-gravitating collapse. 
The Larson-Penston density threshold we adopt for sink creation in our method is a factor 
14 larger than the "Truelove" density threshold adopted in other sink implementations; our 
method also differs from other implementations in that the sink region surrounding each par- 
ticle consists of ghost zones rather than active zones. Outside the sink region, the density 
is below the Truelove value. Our tests show similar results whether we use higher or lower 
density thresholds for sink formation; additional testing in the future can explore whether 
artificial fragmentation is avoided in all cases. We validate our method and implementation 
with a series of tests. These tests include comparison of accretion rates with analytic solu- 
tions for self-similar collapse of isothermal spheres, and comparison of accretion rates and 
solution profiles with ID spherically-symmetric collapse of Bonnor-Ebert spheres and cores 
formed by converging supersonic flows. We demonstrate Galilean invariance of our accretion 
solutions onto sink particles. 

To demonstrate application of our method, we present sample results for a simulation 
of planar converging supersonic flows with turbulence. Filaments form, and local collapse 
produces sink particles which then accrete material from filaments. The mass smoothly 
increases to exceed the mass of the initial bound core in which a sink formed. Sinks forming 
at early stages merge with other smaller sinks if sink regions (of standard size 3 3 zones) 
overlap. In the context of these converging flow tests, we have investigated various criteria 
that have been adopted for creation of sink particles. We find that a (large enough) density 
threshold and local gravitational potential minimum check give the same results as more 
strict sets of criteria in this case (and in other tests we have conducted). 

The implementation of sink particles we present enables study of a wide range of as- 
trophysical applications involving gravitational collapse of a gaseous medium. In particular, 
robust and accurate methods for implementing sink particles make it possible to track se- 
quential formation of multiple protostars over long periods under realistic environmental 
conditions. As such, sink particles represent a crucial numerical tool for addressing key 
unsolved problems in star formation, including the origin of the stellar initial mass func- 
tion. The implementation of sink particles we describe here will be made available to the 
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community in an upcoming release of the Athena code. 
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APPENDIX: VACUUM BOUNDARY CONDITION POTENTIAL VIA 

FOURIER TRANSFORMS 

In this section, we provide details for the solution of Poisson's equation with vacuum 
boundary conditions. Let —4nGQ (x, x') be the Green's function solution for the Poisson 
equation V 2 $ = — 47rG(5(x — x'), so that the potential produced by a density field p(x) is: 

$(x) = -4ttG J £(x,x>(x')d 3 x. (1) 

For potential $ at (x a , y b , z c ) within a domain (L x , L y , L z ) with dimensions (N x , N y , N z ), the 
corresponding discrete sum is: 

N x -lN v -l N g -1 

${x a ,y b , z c ) = -AnG ^2 ^2 S G( x a~ xi,Vb- Vm,z c - z n )p(xi,y m ,z n ); (2) 

1=0 m=0 n=0 

here (a, b, c) and (/, m, n) are the integer indices for the corresponding coordinates x, y, z 
respectively. $(x a , z c ) is a convolution of (?(x) and p(x). 

Because {?(|x — x'|) is a function in the domain [— L X ,L X ] x [— L y ,L y ] x [—L Z ,L Z ], and 
p(x') is a function in the domain [0, L x \ x [0, L y ) x [0, L z ], if we define p(xi,y m , z n ) = for 
/ < 0, m < 0, or n < 0, Equation can be re-written as: 

N x -1 N v -1 N z -1 

®(x a ,y b ,z c ) = -4nG ^ ^ Q{x a - x h y b - y m , z c - z n )p{x h y m , z n ), (3) 

l=-N x m=-N y n=-N z 

We may define periodic functions with period 2L x ,2L y and 2L Z in x, y and z di- 
rections respectively, such that p(x') and £?(x, x') agree with these periodic functions for 
x G [—L x , L x ],y G [—L y , L y ] and z G [—L x , L z }. Then from the Fourier convolution theorem, 
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Equation fl3]) can be expressed in terms of the respective transforms Q and p of Q (x, x') and 
p(x') from the Fourier convolution theorem: 

_4 Q 2N x -l2Ny-\2N z -l fc . 

'^^^' WWW S g E ft^-M* + *Ar), (4) 

Note the summation index for % can be either from — N x — >■ iV^ or — )■ 2A^ — 1 because 
(/ and p are periodic. The same applies to j and k indices. In our implementation within 
Athena, Equation P| is further decomposed into a sum over even and odd terms to save 
memory. 
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Fig. 1. — Circular orbits of two equal mass sink particles orbiting the center of mass at 
different distances, testing our leapfrog particle integrator and TSC particle-mesh Poisson 
solver. The top panels show orbits of diameter d/L = 0.2 at box resolutions (L/Ax) 3 = 
32 3 , 64 3 and 128 3 zones from left to right. The bottom panels show orbits of diameter 
d/L = 0.3 at the same resolutions. Each panel shows the trajectory of one particle for 10 
orbits. 
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Fig. 2. — Accretion history of the centra l sink parti cle for comparison with the case A = 
2.0004 (crosses) and 4 (open circles) from IShul (119771 ). The analytic accretion solutions for 
self-similar flows are shown by the solid lines. The accretion rates from our simulations are 
consistent with the analytical solutions during early stages. The enhancement of accretion 
rates in the later stages is due to the collapse of the outer part of the initial sphere, which is 
affected by boundary conditions. In the units given, the initial sphere masses for the cases 
A = 2.0004 and A = 4 are 18.9 and 40.7, respectively. We fit the linear part of this and 
other accretion histories to derive the accretion rates plotted in Figure |3j 
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Fig. 3. — Accretion rates for self-similar collapse of isothermal spheres with different overden- 
sity coeffici ents, A (se e Equation ({TBI) ). The solid line shows the analytical accretion rates 
derived by IShul ( 119771 ). via direct integration of isothermal fluid equations for self-similar 
models. The solid dots are from accretion rates measured in our 3D simulations. 
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Fig. 4. — Density and velocity field cross-section snapshot for the collapsing near-singular 
isothermal sphere (A = 2.0004). The color scale represents an x — y slice through the density 
(logp). The direction and length of arrows indicate the direction and magnitude of the 
local velocity, with scale as indicated in the upper left. The "expansion wave" is evident in 
the plot, with collapse in the inside (r < 0.67Lj) and a near-static solution in the outside 
(r > 0.67Lj). The location of the expansion wave for the initial conditions was r = 0.43Lj. 
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Fig. 5. — Radial profiles of density (left panel, filled circles) and velocity (right panel, filled 
circles) for evolution of near-singular (A = 2.0004) isothermal sphere. The time interval 
between profiles is constant, equal to 0.l(ir /Gpo) 1 / 2 , with successive profiles moving to lower 
density and velocity over time. In both panel s, th e solid lines are the a nalyt ic solutions 
obtained by integrating Equations (11)-(12) in IShul (119771 ). For the IShul (119771 ) expansion 
wave solution, the outer part of each profile is static (the equilibrium singular isothermal 
sphere), and the inner part of each profile approaches free fall, with p oc r -3 / 2 and v oc r -1 / 2 . 
The dot dash line in the left panel shows the singular isothermal sphere density profile. The 
velocity profiles are plotted in linear-log scale to show the propagation of the "expansion 
wave" . The dashed lines in the right panel indicate the posit i on of the "expansion wave" 
front. The consistency of our numerical solution with the I Shu ( 1977 ) solution is evident. 
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Fig. 6. — Tests of accretion onto a sink particle for the collapsing near-singular isothermal 
sphere (A = 2.0004) with bulk motion across the grid, demonstrating Galilean invariance. 
Horizontal curves, left scale: the temporal evolution of the speed of the sink particle for 
different cases. From bottom to top, the bulk advection speed is 0, 0.5c s , 1.5c s and 2.5c s . 
The time- averaged particle speeds remain constant, and equal to that of the bulk flow. 
Diagonal curves, right scale: the temporal evolution of the mass of the central sink particle 
for the four different cases. The accretion rate is same for static, subsonic, and supersonic 
advection cases, confirming Galilean invariance of our algorithms. 
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Fig. 7. — Density and velocity field cross-section at the instant of singularity for collapse of 
a Bonnor-Ebert sphere. The color scale represents an x — y slice through the volume density 
(logp/p e , for p e the initial density at the outer edge of the sphere). The direction and length 
of arrows indicate the direction and magnitude of the local velocity, with scale as indicated 
in the upper left. The inner part of the sphere is collapsing, with velocity increasing inward 
to approach — 3.4c s . Because we do not have a hot, high-pressure confining medium, the 
outer part of the sphere also expands. 
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Fig. 8. — Comparison of accretion for 3D Cartesian simulation with ID spherically-symmetric 
simulation, for evolution of Bonnor-Ebert sphere. Thick curves, left scale: the temporal 
evolution of the mass accretion rate at the inner edge of the grid. Thin curves, right scale: 
evolution of the central point mass. The solid curves are for the 3D simulation with Athena 
and the dashed curve is for the ID simulation with ZEUS. 



- 34 - 



10000.0 


■ ■ ■ ■ 1 1 1 — ■ — 1 1 ■ — 

: lp 

□s. 


■ ■ ■ ■ i ■ — ■ — ■ — i ■ — -w — ; 


1000.0 






100.0 


- 




10.0 




u u**0^ : 


1.0 






0.1 


■ ■ ■ ■ ' ' 


■ ■ ■ ■ ' ' 



1 



o 

-'> 

-2 
-3 



0.01 



0.10 0.01 0.10 
r/[c s (n/Gp e y/ 2 ] 



Fig. 9. — Density (left) and radial velocity (right) profile comparison between 3D Athena 
simulation and ID spherically-symmetric ZEUS simulation for the instant of singularity in 
Larson-Penston collapse. The pluses are for the ID simulation and the squares are for the 
3D simulation. Also shown is the singular LP density profile (Equation ([7j); solid line on 
left). Units are normalized using the initial density at the edge of the sphere, p e . The initial 
outer radius of the sphere is l.lrBE.crit, and the initial density in the region exterior to the 
sphere is 0.001p e . 
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Fig. 10. — Radial density (left) and velocity (right) profiles throughout collapse and infall 
stages for test beginning with static Bonnor-Ebert density profile. The solid lines are during 
collapse (i = for the lower and t = 0.043to for the upper curve in the left panel, and the 
opposite in the right panel), the dashed lines are at the instant of sink particle creation, and 
the dotted lines are during the infall stage (t = 0.129t for the upper and t = 0.172t for 
the lower curve in each panel). The time unit is t = c s (7i/Gp e Y^ 2 for p e the initial density 
at the edge of the sphere. The behavi or in each stage is cons istent with previous results 
from spherically-symmetric simulations (IGong fc Ostrikerl 120091 ) . The time interval between 
profiles is At = 0.043t - 
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Fig. 11. — Comparison of accretion evolution for 3D Cartesian and for ID spherically- 
symmetric simulation for spherical converging supersonic (Mach 2) flow. Heavy curves show 
accretion rate, and light curves show sink particle mass. Solid curves are the 3D model with 
Athena, and dashed curves are the ID model with ZEUS. 
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Fig. 12. — Evolution of surface density projected in the ^-direction (color scale log£/E ) 
for a planar converging flow simulation with inflow Mach number hA = 5, and supersonic 
turbulence. The four panels from top left to bottom right show snapshots at four instants: 
0.301tj, 0.349tj, 0.398tj, and 0.446tj. The top left panel shows the surface density when the 
first sink "1" is created. The numbers in panels mark the time sequence of sink formation. 
The black solid dots in the first three frames show instantaneous locations of sink particles. 
In the final panel, the four white large solid dots show the surviving sinks after mergers. The 
magenta curves shows the trajectories of all sinks from birth (marked with triangles) to the 
end of the simulation. The area of the projected domain is Lj x Lj (see Equation (1201)). 
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Fig. 13. — Temporal evolution of sink particle masses (see Equation (12T1) and Equation (1221) 
for definition of Mj and tj respectively). The solid black lines show the sinks that survive to 
late times. The dotted lines show the sinks that are eventually merged into larger particles. 
The solid black dots show the gravitationally bound dense core masses immediately before 
the creation of the corresponding sink particle at its center. 
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Fig. 14. — Temporal evolution of sink particle masses based on adoption of different creation 
criteria. The solid lines are for sink particles created with the LP density threshold, and 
local potential minimum criteria. The tracks marked by pluses are for sink particles created 
based on additional criteria: a converging flow and a gravitationally bound state. The 
dotted lines are for sink particles created based on the Truelove density threshold and the 
local gravitational potential minimum check. 




